program xyz2para;
uses common,decor,ppot,tools,inp,outp;
Type
   r_xyz  = array[1..nsites] of v3;
   r_e	  = array[1..nsites] of myreal;
   r_tp	  = array[1..nsites] of integer;
Var
   i,j,k,l,m,n,o,dbg,j1,j2    : integer;
   txt			      : text;
   r,rl			      : r_xyz;
   ip			      : array[1..nsites] of atomrec;
   rb			      : array[1..nsites] of v33;
   ee			      : r_e;
   itp,ito		      : r_tp;
   nsteps,omax,tinp,natr,nats : integer;
   ich			      : array[1..atps] of integer;
   xgamma,dr2,drmax	      : myreal;
   duo,duos		      : array[1..t_orb,1..atps] of v2;
   no			      : array[1..t_orb,1..atps] of integer;
   tname,xyzname	      : string[50];
   rbt			      : v33;
   w1,w2,w3,rx		      : v3;
   nofr			      : integer;
   ofr			      : array[1..t_orb] of integer;
   nnele		      : integer;
   ele			      : i10;
   dr_max		      : myreal;
   po_max		      : i3;

procedure dparinit;
var i,j,o : integer;
begin
   for i:=1 to n_pxyz do begin
      for j:=1 to nn_pxyz[i] do begin ndps[i,j]:=0; d_pxyz[i,j]:=0.0; end;
   end;
end;
procedure dpara2xyz(rec : atomrec;ch:integer; Var x:v3);
{NOTICE: using globals nn_pxyz, ch_pxyz, dpxyz}
Var i,j,o,po,ma,ma0,pox : integer;
   done			 : boolean;
begin
   for o:=1 to 3 do x[o]:=0.0;
   {check we have chemical parameter for given atom}
   ma:=0; ma0:=0; for i:=1 to 3 do if rec[i]>0 then begin
      ma0:=ma0+1;
      for j:=1 to nn_pxyz[rec[i]] do if ch_pxyz[rec[i],j]=ch then ma:=ma+1;
   end;
   {NOTE: case ma0=0 for completely fixed orbit is treated automatically}
   {add new chemistry}
   if (ma0<>ma) then begin
      if ma>0 then begin
	 writeln('Atom-chemistry parameter error: chemistry cannot be different for different cartezian components!');
	 for o:=1 to 6 do write(rec[i]:1,' '); writeln;
	 halt;
      end;
      add_chem_para(rec,ch);
   end;
   for i:=1 to 3 do if rec[i]>0 then begin
      po:=rec[i]; {parameter}
      pox:=rec[i+3]; {basis pointer}
      if i<3 then begin {deca-plane}
	 {chemistry-specific displacement}
	 for j:=1 to nn_pxyz[po] do if (ch=ch_pxyz[po,j]) then begin
	    for o:=1 to 2 do x[o]:=x[o]+d_pxyz[po,j]*re[pox,o];
	 end;
      end else begin {periodic direction}
	 {chemistry-specific displacement}
	 for j:=1 to nn_pxyz[po] do if (ch=ch_pxyz[po,j]) then begin
	    x[i]:=x[i]+d_pxyz[po,j]*pox;
	 end;
      end;
   end;
end;
procedure xyz2pars(nat : integer; itp:r_tp; Var r,rl:r_xyz; Var dr2 : myreal);
Var
   i,j,k,o,m,n	      : integer;
   bbij,ar2,cf	      : myreal;
   drt,drc	      : v3;
   bb		      : array[1..nsites] of v33;
   io,jo	      : integer;
   i1,i2,p1,p2,p3,ina : integer;
   p		      : i3;
   ar,u,ud,udd,drm2   : myreal;
   v,w,w1,drx	      : v3;
   bo		      : boolean;
   b,dra	      : r_xyz;
   z1,z2	      : v2;
   update	      : boolean;
begin
   drm2:=drmax*drmax;
   epot:=0.0; ina:=0; dr2:=0.0;
   dparinit;
   dr_max:=0.0;
   for o:=1 to 3 do po_max[o]:=0;
   for i:=1 to nat do begin
      for o:=1 to 3 do drc[o]:=bound(rl[i,o]-r[i,o]); {dr-cryst}
      for m:=1 to 3 do drt[m]:=0.0;
      for m:=1 to 3 do
	 for n:=1 to 3 do drt[n]:=drt[n]+drc[m]*bl3[m,n]; {cryst2cart}
      ar2:=0.0; for o:=1 to 3 do ar2:=ar2+drt[o]*drt[o];
      if (ar2>dr_max) then begin
	 {writeln(i:1,' ',sqrt(ar2):1:5);}
	 dr_max:=sqrt(ar2);
	 po_max[1]:=ito[i];
      end;
      dr2:=dr2+ar2;
      mxle3(blt,drt,dra[i]);
      mxle3(rb[i],drt,w1);
      if (ofr[ito[i]]>0) then add2para(ip[i],itp[i],w1);
   end;
   {writeln(n_pxyz:1);}
   {update parameters}
   for i:=1 to n_pxyz do
      for j:=1 to nn_pxyz[i] do if (ndps[i,j]>0) then begin
	 d_pxyz[i,j]:=d_pxyz[i,j]/ndps[i,j];
	 {writeln(i:1,' ',j:1,' ',ndps[i,j]:1,' ',dpxyz[i,j]:1:5,' ',d_pxyz[i,j]:1:5);}
	 dpxyz[i,j]:=dpxyz[i,j]+d_pxyz[i,j];
      end;
   for i:=1 to nat do begin {update coordinates}
      dpara2xyz(ip[i],itp[i],w);
      for o:=1 to 2 do z1[o]:=w[o];
      tra(z1,z2);
      for o:=1 to 2 do drc[o]:=z2[o];
      drc[3]:=bound(w[3]/multi); {now drc is SYMMETRIC du in cryst format}
      {orbit's mean SYMMETRIC displacement}
      for m:=1 to 3 do drt[m]:=0.0;
      for m:=1 to 3 do
	 for n:=1 to 3 do drt[n]:=drt[n]+drc[m]*bl3[m,n]; {cryst2cart}
      ar2:=0.0; for o:=1 to 2 do ar2:=ar2+drt[o]*drt[o];
      duos[ito[i],itp[i],1]:=duos[ito[i],itp[i],1]+ar2; {deca}
      duos[ito[i],itp[i],2]:=duos[ito[i],itp[i],2]+drt[3]*drt[3]; {per}
      for o:=1 to 3 do r[i,o]:=bound(r[i,o]+drc[o]); {symmetric cryst R}
      {orbit's mean SYMMETRY-BREAKING displacement}
      ar2:=0.0;
      for o:=1 to 3 do drc[o]:=bound(rl[i,o]-r[i,o]); {cryst}
      for m:=1 to 3 do drt[m]:=0.0;
      for m:=1 to 3 do
	 for n:=1 to 3 do drt[n]:=drt[n]+drc[m]*bl3[m,n]; {cryst2cart}
      ar2:=0.0; for o:=1 to 2 do ar2:=ar2+drt[o]*drt[o];
      duo[ito[i],itp[i],1]:=duo[ito[i],itp[i],1]+ar2; {deca}
      duo[ito[i],itp[i],2]:=duo[ito[i],itp[i],2]+drt[3]*drt[3]; {per}
      no[ito[i],itp[i]]:=no[ito[i],itp[i]]+1; {count}
   end;
end; { xyz2pars }
begin
   for o:=1 to t_orb do ofr[o]:=1;
   {need before read_deco can be read!!}
   reset(txt,'INPUT/objects.def');
   readln(txt,deconame);
   readln(txt,nu_offs);
   readln(txt,dpmodel);
   readln(txt,o); repre5:=(o=1);
   close(txt);
   readln(tinp);
   if(tinp>1)then begin
      writeln('tinp=',tinp:1,' : not supported');
      exit;
   end;
   readln(tname);
   readln(xyzname);
   dbg:=0;
   BASIS;
   reset(txt,tname);
   IMPORT_TILING_BODY(tinp,txt);
   TABLE_INIT;
   if (tinp>0) then IMPORT_CONFIG(txt);
   close(txt);
   UNIT_CELL(multi);  
   for i:=1 to 3 do begin
      for o:=1 to 3 do blt[o,i]:=bl3[i,o];
   end;
   PP_INP;
   if (tinp=0) then begin
      i:=0;
      READ_DECO(i); {reads objects.def, atoms.def}
      for i:=1 to tnodes do nan[i]:=0;
      for i:=1 to tnodes do NEW_NODE_OBJECTS(i,1);
   end;
   omax:=0;
   nat:=0;
   {writeln('DBG sigma',sigmaflip);}
   write('# Initializing atomic decoration from '+tname,'(input_type=',tinp:1,'): ');
   for i:=1 to tnodes do begin
      for j:=1 to nan[i] do begin nat:=nat+1;
	 for o:=1 to 3 do r[nat,o]:=bound(n_r[i,j,o]+rt[i,o]);
	 itp[nat]:=n_a[i,j,1];
	 ito[nat]:=n_a[i,j,2];
	 ip[nat]:=n_para[i,j];
	 add_chem_para(n_para[i,j],itp[nat]);
	 {basis for each atom relative to object origin [node]}
	 get_atom_basis(ip[nat],rb[nat]);
	 if (ito[nat]>omax) then omax:=ito[nat];
	 if (itp[nat]=0) then nat:=nat-1 else begin
	    {write(nat:1,' ');
	    for o:=1 to 3 do write(r[nat,o]:1:5,' '); writeln('  ',itp[nat]:1);}
	 end;
      end;
   end;
   writeln(nat:1,' atoms.');
   writeln('# Reading xyz_rel file : ',xyzname);
   reset(txt,xyzname);
   for o:=1 to 3 do readln(txt); {basis}
   readln(txt,nats);
   natr:=0;
   for i:=1 to nats do begin
      natr:=natr+1;
      for o:=1 to 3 do read(txt,rx[o]); readln(txt,j1,j2);
      if (j1=0) then natr:=natr-1 else begin
	 {if not(j1=itp[natr]) or not(j2=ito[natr]) then begin}
	 if not(j1=itp[natr]) then begin
	    writeln('Incompatible xyz/config: ',j1:1,' <> ',itp[natr]:1);
	    halt;
	 end;
	 for o:=1 to 3 do rl[natr,o]:=rx[o];
      end;
   end;
   if (natr<>nat) then begin
      writeln('nat<>natr : ',nat:1,' ',natr:1,'  HALTING');
      halt;
   end;
   close(txt);
   for i:=1 to atps do ich[i]:=0;
   for i:=1 to nat do ich[itp[i]]:=ich[itp[i]]+1;
   for i:=1 to t_orb do for j:=1 to atps do begin  {orbit energies}
      for o:=1 to 2 do begin
	 duo[i,j,o]:=0.0; duos[i,j,o]:=0.0;
      end;
      no[i,j]:=0;
   end;   
   xyz2pars(nat,itp,r,rl,dr2);

   rewrite(txt,'xyz.apa');
   for i:=1 to 3 do begin
      for j:=1 to 3 do write(txt,bl3[i,j]:12:5); writeln(txt);
   end;
   writeln(txt,nat:1,' # atoms (no vacancies)');
   for i:=1 to nat do begin
      for o:=1 to 3 do write(txt,r[i,o]:15:8);
      writeln(txt,itp[i]:6,ito[i]:6);
   end;
   close(txt);
   export_config('copa.apa');
   writeln('### ORB CHEM   dR_sym dZ_sym  dR_break dZ_break [A]');
   for i:=1 to omax do for j:=1 to atps do if (no[i,j]>0) then begin
      write(i:2,' ',j:2,'    ');
      write(sqrt(duos[i,j,1]/no[i,j]):1:4,'  ',sqrt(duos[i,j,2]/no[i,j]):1:4);
      writeln('   ',sqrt(duo[i,j,1]/no[i,j]):1:4,'  ',sqrt(duo[i,j,2]/no[i,j]):1:4);
   end;
end.

